Setup MKL direct factorizations #349

ChrisRackauckas · 2023-07-31T13:29:10Z

MWE:

using LinearSolve, MKL_jll
A = rand(4, 4); b = rand(4); u0 = zeros(4);
lp = LinearProblem(A, b);
truesol = solve(lp, LUFactorization())
mklsol = solve(lp, MKLLUFactorization())
@test truesol ≈ mklsol

The segfault can be reproduced just with the triangular solver. MWE without LinearSolve:

using MKL_jll
using LinearAlgebra: BlasInt, LU
using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, chkargsok
const usemkl = MKL_jll.is_available()

function getrf!(A::AbstractMatrix{<:Float64}; ipiv = similar(A, BlasInt, min(size(A,1),size(A,2))), info = Ref{BlasInt}(), check = false)
    require_one_based_indexing(A)
    check && chkfinite(A)
    chkstride1(A)
    m, n = size(A)
    lda  = max(1,stride(A, 2))
    ccall((:dgetrf_, MKL_jll.libmkl_rt), Cvoid,
            (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64},
            Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
            m, n, A, lda, ipiv, info)
    chkargsok(info[])
    A, ipiv, info[] #Error code is stored in LU factorization type
end

A = rand(4,4); b = rand(4)
getrf!(A)
LU(getrf!(A)...) \ b

@test

MWE: ```julia using LinearSolve, MKL_jll A = rand(4, 4); b = rand(4); u0 = zeros(4); lp = LinearProblem(A, b); truesol = solve(lp, LUFactorization()) mklsol = solve(lp, MKLLUFactorization()) @test truesol ≈ mklsol ``` The segfault can be reproduced just with the triangular solver. MWE without LinearSolve: ```julia using MKL_jll using LinearAlgebra: BlasInt, LU using LinearAlgebra.LAPACK: require_one_based_indexing, chkfinite, chkstride1, chkargsok const usemkl = MKL_jll.is_available() function getrf!(A::AbstractMatrix{<:Float64}; ipiv = similar(A, BlasInt, min(size(A,1),size(A,2))), info = Ref{BlasInt}(), check = false) require_one_based_indexing(A) check && chkfinite(A) chkstride1(A) m, n = size(A) lda = max(1,stride(A, 2)) ccall((:dgetrf_, MKL_jll.libmkl_rt), Cvoid, (Ref{BlasInt}, Ref{BlasInt}, Ptr{Float64}, Ref{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), m, n, A, lda, ipiv, info) chkargsok(info[]) A, ipiv, info[] #Error code is stored in LU factorization type end A = rand(4,4); b = rand(4) getrf!(A) LU(getrf!(A)...) \ b ```

ai-maintainer

AI-Maintainer Review for PR - WIP: Setup MKL direct factorizations

Title and Description 👍

The Title and Description are clear but could be more detailed

The title and description of the pull request do provide some explanation of the purpose of the changes. The title indicates that it is a work in progress (WIP) for setting up MKL direct factorizations. The description includes a minimal working example (MWE) that demonstrates the issue and provides code snippets to reproduce the segfault. However, the description could be more detailed and provide additional context about why these changes are necessary and what problem they aim to solve.

Scope of Changes 👍

The changes are narrowly focused

The changes in this pull request appear to be narrowly focused on setting up MKL direct factorizations. The author is not trying to resolve multiple issues simultaneously. The changes are specific to integrating MKL functionality into the LinearSolve package and addressing the segfault that occurs with the triangular solver.

Documentation 👍

All new functions, classes, or methods have docstrings

Yes, all newly added functions, classes, or methods have docstrings describing their behavior, arguments, and return values. There are no missing docstrings in this pull request.

Testing 👍

The author has provided a test case, but more details would be helpful

Yes, the description includes a description of how the author tested the changes. It provides a minimal working example (MWE) that demonstrates the issue and includes code snippets to reproduce the segfault. The author also includes a test assertion (`@test truesol ≈ mklsol`) to verify that the MKL-based solution is approximately equal to the true solution. However, it would be helpful if the author provided more details about the testing process, such as any additional test cases or scenarios that were considered.

Suggested Changes

Please provide more context in the PR description about why these changes are necessary and what problem they aim to solve.
Please provide more details about the testing process, such as any additional test cases or scenarios that were considered.

Reviewed with AI Maintainer

ChrisRackauckas · 2023-07-31T13:29:34Z

@staticfloat or @chriselrod maybe you might have an idea of what I'm doing wrong here to cause a segfault?

codecov · 2023-07-31T13:33:46Z

Codecov Report

Merging #349 (aaf64d3) into main (98a2292) will decrease coverage by 10.60%.
The diff coverage is 88.46%.

@@             Coverage Diff             @@
##             main     #349       +/-   ##
===========================================
- Coverage   75.29%   64.70%   -10.60%     
===========================================
  Files          17       18        +1     
  Lines        1267     1292       +25     
===========================================
- Hits          954      836      -118     
- Misses        313      456      +143

Files Changed	Coverage Δ
src/LinearSolve.jl	`84.09% <50.00%> (-1.63%)`	⬇️
ext/LinearSolveMKLExt.jl	`91.30% <91.30%> (ø)`
src/extension_algs.jl	`12.50% <100.00%> (-57.07%)`	⬇️

... and 4 files with indirect coverage changes

📣 We’re building smart automated test selection to slash your CI/CD build times. Learn more

chriselrod · 2023-07-31T15:16:05Z

ext/LinearSolveMKLExt.jl

+        ipiv = similar(A, BlasInt, min(size(A,1),size(A,2)))
+    end
+
+    ccall((:dgetrf_, MKL_jll.libmkl_rt), Cvoid,


https://github.com/JuliaLinearAlgebra/BLASBenchmarksCPU.jl/blob/8d26d05539d675ed91f12ce7a14007b87ed01e51/src/ccallblas.jl#L71

I don't seem to be adding a trailing underscore to dgetrf_

staticfloat · 2023-07-31T23:56:51Z

You've got MKL in 32-bit mode, but passing it Int64's, so the pivot indices that come back are gibberish, which causes the future LU factorization to have bogus indices, and the \ to segfault:

julia> A, ipiv, info = getrf!(A);

julia> ipiv
4-element Vector{Int64}:
  8589934595
 17179869187
           0
           0

There are a bunch of ways to fix this, but the best is to import @blasfunc (remember that guy we dropped when copy-pasting) from LinearAlgebra.BLAS and use that via @blasfunc(dgetrf_). That's good because you're already polymorphic by using similar(..., BlasInt, ...) when initializing ipiv, so this just makes it match up to the correct version by default.

ai-maintainer bot reviewed Jul 31, 2023

View reviewed changes

chriselrod reviewed Jul 31, 2023

View reviewed changes

fix and test

16c775e

ChrisRackauckas changed the title ~~WIP: Setup MKL direct factorizations~~ Setup MKL direct factorizations Aug 1, 2023

ChrisRackauckas added 3 commits July 31, 2023 23:47

fix MKL test

d00d025

add missing using

3c1e66d

add backwards compat

aaf64d3

ChrisRackauckas merged commit 1ee467b into main Aug 1, 2023
14 of 17 checks passed

ChrisRackauckas deleted the mkl branch August 1, 2023 05:56

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Setup MKL direct factorizations #349

Setup MKL direct factorizations #349

ChrisRackauckas commented Jul 31, 2023

ai-maintainer bot left a comment •

edited

Loading

ChrisRackauckas commented Jul 31, 2023

codecov bot commented Jul 31, 2023 •

edited

Loading

chriselrod Jul 31, 2023 •

edited

Loading

staticfloat commented Jul 31, 2023 •

edited

Loading

Setup MKL direct factorizations #349

Setup MKL direct factorizations #349

Conversation

ChrisRackauckas commented Jul 31, 2023

ai-maintainer bot left a comment • edited Loading

Choose a reason for hiding this comment

AI-Maintainer Review for PR - WIP: Setup MKL direct factorizations

Title and Description 👍

Scope of Changes 👍

Documentation 👍

Testing 👍

Suggested Changes

ChrisRackauckas commented Jul 31, 2023

codecov bot commented Jul 31, 2023 • edited Loading

Codecov Report

chriselrod Jul 31, 2023 • edited Loading

Choose a reason for hiding this comment

staticfloat commented Jul 31, 2023 • edited Loading

ai-maintainer bot left a comment •

edited

Loading

codecov bot commented Jul 31, 2023 •

edited

Loading

chriselrod Jul 31, 2023 •

edited

Loading

staticfloat commented Jul 31, 2023 •

edited

Loading